NASA/TP-2007-215086 



Refinement of Timoshenko Beam Theory for 
Composite and Sandwich Beams Using Zigzag 
Kinematics 


Alexander Tessler 

Langley Research Center, Hampton, Virginia 

Marco Di Sciuva and Marco Gherlone 
Department of Aeronautics and Space Engineering 
Politecnico di Torino, Torino, Italy 


November 2007 



The NASA STI Program Office ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA Scientific and Technical Infonnation (STI) 
Program Office plays a key part in helping NASA 
maintain this important role. 

The NASA STI Program Office is operated by 
Langley Research Center, the lead center for NASA’s 
scientific and technical infonnation. The NASA STI 
Program Office provides access to the NASA STI 
Database, the largest collection of aeronautical and 
space science STI in the world. The Program Office is 
also NASA’s institutional mechanism for 
disseminating the results of its research and 
development activities. These results are published by 
NASA in the NASA STI Report Series, which 
includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
programs and include extensive data or 
theoretical analysis. Includes compilations of 
significant scientific and technical data and 
information deemed to be of continuing 
reference value. NASA counterpart of peer- 
reviewed fonnal professional papers, but having 
less stringent limitations on manuscript length 
and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or of 
specialized interest, e.g., quick release reports, 
working papers, and bibliographies that contain 
minimal annotation. Does not contain extensive 
analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical infonnation from NASA 
programs, projects, and missions, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific and 
technical material pertinent to NASA’s mission. 

Specialized services that complement the STI 
Program Office’s diverse offerings include creating 
custom thesauri, building customized databases, 
organizing and publishing research results ... even 
providing videos. 

For more infonnation about the NASA STI Program 
Office, see the following: 

• Access the NASA STI Program Home Page at 
http://www. sti. nasa.gov 

• E-mail your question via the Internet to 
help@sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 
at (301) 621-0134 

• Phone the NASA STI Help Desk at 
(301)621-0390 

• Write to: 

NASA STI Help Desk 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 



NASA/TP-2007-215086 



Refinement of Timoshenko Beam Theory for 
Composite and Sandwich Beams Using Zigzag 
Kinematics 


Alexander Tessler 

Langley Research Center, Hampton, Virginia 

Marco Di Sciuva and Marco Gherlone 
Department of Aeronautics and Space Engineering 
Politecnico di Torino, Torino, Italy 


National Aeronautics and 
Space Administration 

Langley Research Center 
Hampton, Virginia 23681-2199 


November 2007 



The use of trademarks or names of manufacturers in this report is for accurate reporting and does not 
constitute an official endorsement, either expressed or implied, of such products or manufacturers by the 
National Aeronautics and Space Administration. 


Available from: 


NASA Center for AeroSpace Information (CASI) 
7115 Standard Drive 
Hanover, MD 21076-1320 
(301) 621-0390 


National Technical Information Service (NTIS) 
5285 Port Royal Road 
Springfield, VA 22161-2171 
(703) 605-6000 



Abstract 


A new refined theoiy for laminated-composite and sandwich beams 
that contains the kinematics of the Timoshenko Beam Theoiy as a proper 
baseline subset is presented. This variationally consistent theory is 
derived from the virtual work principle and employs a novel piecewise 
linear zigzag function that provides a more realistic representation of the 
deformation states of transverse shear flexible beams than other similar 
theories. This new zigzag function is unique in that it vanishes at the top 
and bottom bounding surfaces of a beam. The formulation does not 
enforce continuity of the transverse shear stress across the beam ’s cross- 
section, yet is robust. Two major shortcomings that are inherent in the 
previous zigzag theories, shear-force inconsistency and difficulties in 
simulating clamped boundary conditions, and that have greatly limited 
the utility of these previous theories are discussed in detail. An approach 
that has successfully resolved these shortcomings is presented herein. 
Exact solutions for simply supported and cantilevered beams subjected 
to static loads are derived and the improved modelling capability of the 
new “zigzag” beam theory is demonstrated. In particular, extensive 
results for thick beams with highly heterogeneous material lay-ups are 
discussed and compared with corresponding results obtained from 
elasticity solutions, two other “zigzag” theories, and high-fidelity finite 
element analyses. Comparisons with the baseline Timoshenko Beam 
Theoiy are also presented. The comparisons clearly show the improved 
accuracy of the new, refined “zigzag” theoiy presented herein over 
similar existing theories. This new theory can be readily extended to 
plate and shell structures, and should be useful for obtaining relatively 
low-cost, accurate estimates of structural response needed to design an 
important class of high-performance aerospace structures. 


Nomenclature 


A, 2 h, L 

beam’s cross-sectional area, depth, and span 

2 h (k) 

thickness of the k -th layer 

q,F 

transverse pressure [force/length] and force loading 


axial and shear tractions prescribed at the ends of beam x (q = a,b ) 

(uf, „!*>) 

axial and transverse components of the displacement vector in the 
k th layer 

(u, <9, 1 //, w) 

kinematic variables of zigzag theory 
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zigzag function 
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axial and transverse shear strain in the k -th layer 
axial and transverse shear stress in the k -th layer 
axial and shear modulus of the k -th layer 
penalty factors 

local ( k -th layer) thickness coordinate 

normalized axial displacement along the interface between 
the k -th and ( k + 1) layers 

stress resultants of the refined zigzag beam theory 

constitutive stiffness coefficients of the refined zigzag beam theory 
number of material layers through the beam’s depth 
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1. Introduction 


Performance and weight advantages of advanced composite materials have led to their sustained and 
increased application to military and civilian aircraft, aerospace vehicles, naval, and civil structures. To 
design efficient and reliable laminated-composite structures, improved analytical and computational 
methods that accurately incorporate the principal non-classical effects are necessary. In relatively thick 
beams, and even in thinner beams that are highly heterogeneous, transverse shear deformation may 
substantially influence design-significant response quantities such as normal stresses, deflections, and 
vibration modes and natural frequencies. For this class of problems, classical deformation theories are 
generally less than adequate. This shortcoming is particularly pronounced in relatively thick laminated- 
composite and sandwich structures with material layers that exhibit large differences in the transverse 
shear properties. Application of classical theories in this context often leads to non-conservative 
predictions for deformation, stresses, and natural frequencies. Moreover, in classical shear deformation 
models, transverse shear stresses fail to satisfy equilibrium conditions at the layer interfaces, resulting in 
excessively erroneous discontinuities in these quantities. The two fundamental assumptions of the 
classical Bemoulli-Euler Beam Theory are that the transverse shear and through-the-depth normal strains 
are negligible, compared to the axial strain associated with bending action. These assumptions have been 
shown to be valid for the bending of very slender beams. As a result, the bending deformation may be 
characterized by a single deflection variable. However, Hooke’s law yields a zero-valued transverse shear 
stress, which violates equilibrium of the internal forces acting within a beam. To rectify this inconsistency 
in the classical Bemoulli-Euler Beam Theory, a beam equilibrium equation is used to obtain the internal 
transverse shear force from which an average shear stress is computed. Timoshe nk o [1] derived a new 
beam theory by adding an additional kinematic variable in the displacement assumptions, the bending 
rotation, and a shear correction factor to provide a first-order account of shear deformation in an average 
sense, while retaining the assumption that the through-the-depth normal strain is negligible. This 
modification to the kinematics allows the transverse shear stresses to be obtained from Hooke’s law and 
extends the range of applicability of the resulting theory to moderately thick beams and higher vibration 
frequencies associated with short wave length modes. 

The Timoshenko Beam Theory (TBT) and analogous shear-deformation theories for plate and shell 
structures have been widely used in structural analysis of homogeneous and laminated-composite 
structures. This theory produces inadequate predictions, however, when applied to relatively thick 
laminated composites with material layers that have highly dissimilar stiffness characteristics. Even with 
a judiciously chosen shear correction factor, which is dependent on the laminate stacking sequence, the 
TBT tends to underestimate, often substantially, the axial stresses at the top and bottom, outer fibers of a 
beam. Moreover, along the layer interfaces of a laminated-composite beam, the transverse shear stresses 
predicted often exhibit excessively erroneous discontinuities. The reason for these difficulties is likely 
associated with the higher complexity of the “true” variation of the displacement field across a highly 
heterogeneous beam cross-section. Clearly, the linear through-thickness displacement assumption for the 
axial displacement is the main shortcoming of the TBT when it is applied to complex material systems 
with significant transverse shear flexibility. 

The need for beam, plate, and shell theories with better predictive capabilities has led to the 
development of the so-called higher-order theories that are also commonly known as refined equivalent 
single-layer theories (for a review of many of such theories, refer to [2]). In these refined theories, 
higher-order kinematic terms, with respect to the beam-depth or shell-thickness coordinate, have been 
added to the expressions for the in-plane displacements and, in some cases, to the expressions for the 
transverse displacement. While notable response improvements have been achieved with several of these 
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theories, they are generally inadequate for predicting the correct shear and axial stresses in, and the high- 
frequency dynamic responses of, highly heterogeneous laminated-composite and sandwich structures that 
are moderately thick. Because of these shortcomings, refined theories that incorporate deformation 
characteristics of the individual layers of the structure have been developed and are commonly referred to 
as layer-wise theories. Notable early contributions to these theories were made by Ambartsumian [3] and 
by Sun and Whitney [4], 

In general, layer-wise theories presume that the laminate deformation and stress fields depend 
significantly on the behavior of the individual layers. In these theories, the layer kinematics are 
independently described and certain physical continuity requirements are enforced (e.g., refer to the 
review article in [2]). More specifically, the layer-wise representation of the kinematics yields zigzag-like 
displacements through the beam depth or shell thickness. Moreover, the increased kinematic freedom 
enables the enforcement of stress continuity conditions at the layer interfaces. The major drawback of 
these theories, however, is that the number of kinematic variables is dependent on the number of layers in 
the beam or shell. This attribute typically results in a great number of solution variables which makes the 
general layer-wise theories computationally unattractive and, as a result, detracts from their usefulness. 
Moreover, these theories are particularly cumbersome to implement within a displacement-based finite 
element method. For further discussions on the merits and drawbacks of general layer-wise theories, refer 
to [5], 

The so-called zigzag theories are an important special sub-class of the general layer-wise theories that 
appear to be amenable to engineering practice. These theories also assume a zigzag pattern for the in- 
plane displacements and enforce continuity of the transverse shear stresses across the entire laminated- 
beam depth or laminated-shell thickness. Importantly, the number of kinematic variables in these special 
layer-wise theories is independent of the number of layers. Early works on the zigzag theories were 
presented by Di Sciuva [6-8] and Murakami [9] for the bending of plates. Di Sciuva [10-11] also 
demonstrated that zigzag theories are generally well-suited for finite element implementation. 

In Di Sciuva’s earlier efforts on plates [6-7], a piecewise linear inplane-displacement function is added 
to a first-order shear deformation theory (FSDT) in which the transverse shear strains are used as the 
primary kinematic variables instead of the bending rotations. This piecewise linear function, and others 
like it, are referred to herein as a “zigzag” function. To retain only the kinematic variables of the FSDT, 
Di Sciuva enforced a constant shear stress across the entire laminate thickness. This procedure led to the 
desired enhancement in the representation of the inplane displacement fields and simultaneously achieved 
transverse shear stress continuity along layer interfaces. An important characteristic of Di Sciuva’s 
formulation is that, when the zigzag function vanishes for homogeneous plates, the theory is reduced to 
the FSDT baseline. Later, Di Sciuva [12-13] introduced additional enhancements to his zigzag model by 
adding a cubic in-plane displacement to the zigzag function. It is important to note that the Di Sciuva 
theories require (^-continuous shape functions for formulating suitable finite elements. These 
approximation schemes are significantly less attractive, especially for plate and shell finite elements, than 
the (^-continuous displacement interpolation schemes associated with FSDTs. 

With the goal of obtaining ('^-continuous finite elements, Averill explored the new linear [14] and 
cubic [15] zigzag beam models. In particular, Averill modified Di Sciuva’s approach by using the TBT as 
his underlying baseline theory, adding an additional kinematic variable associated with a zigzag function, 
and by introducing an ad hoc penalty term to the variational principle. The penalty term enforces 
continuity of the transverse shear stress across the beam cross-section in a limiting sense as the penalty 
parameter approaches a large value. 

Di Sciuva’s theory runs into theoretical difficulties with regard to the physical significance of the 
transverse shear stress associated with the theory. The difficulty is especially evident at a clamped 
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support, where the cross-sectional area integral of the shear stress, obtained from constitutive relations, 
does not correspond to the total shear force required by equilibrium. Thus, the correct shear force and the 
average shear stress must be determined from the equilibrium equation that relates the shear force to the 
derivative of the bending moment, as in the Bemoulli-Euler Beam Theory. AverilTs theory also suffers 
from its inability to model correctly a clamped boundary condition, where it predicts erroneously that the 
transverse shear stress and the corresponding resultant force vanish. To alleviate this anomaly, Averill 
proposed a boundary condition compromise at the expense of variational consistency of the theory, in 
which a kinematic variable representing the amplitude of the zigzag displacement is left out of the 
variationally required set of boundary conditions. Consequently, extensive analytic and numerical studies 
that have been conducted primarily focused on beams and plates with simply supported boundaries [6, 7, 
and 12-15]. Recently, a zigzag plate analysis was discussed for clamped plates [16]; however, no results 
were presented for the shear stresses at a clamped edge. 

Scrutiny of the zigzag theories discussed herein has revealed some serious shortcomings. The 
objective of the present study is to present a new refined zigzag theory that is free of these shortcomings 
and amenable to finite element implementation. In particular, the present paper discusses a new refined 
zigzag beam theory that is consistently derived from the virtual work principle, by refining the ideas of 
Timoshenko, Di Sciuva, and Averill. The key attributes of the present theory are, first, the proposed 
zigzag function vanishes at the top and bottom surfaces of the beam and does not require full shear-stress 
continuity across the laminated-beam depth. Second, all boundary conditions, including the fully clamped 
condition, can be modeled adequately. And third, the theory only requires C°-continuous kinematics for 
finite element modeling, as do elements based on the theories of Timoshe nk o [1], Mindlin [17], and 
Reissner [18]. This latter attribute lends itself to developing computationally efficient and robust beam, 
plate, and shell finite elements. Overall, the theory appears as a natural extension of the TBT to 
laminated-composite beams, and it is devoid of the drawbacks of the zigzag theories discussed previously. 

In the remainder of the report, the concept of zigzag kinematic assumptions is first described. Then the 
zigzag schemes of Di Sciuva and Averill are elaborated in detail, and their deficiencies with respect to the 
transverse shear properties and clamped boundary conditions are highlighted. A unique zigzag function is 
then introduced to formulate the basis for the refined zigzag theory, giving rise to the transverse shear 
stress that has a piecewise constant distribution across the laminate thickness. As an added explanation of 
the underlying reasons for the drawbacks of AverilTs formulation, a penalized form of the constitutive 
equations is introduced within the present theory. The equations of equilibrium and associated boundary 
conditions are then derived from the virtual work principle. Finally, the refined zigzag theory is assessed 
quantitatively by way of exact solutions for simply supported and cantilevered composite and sandwich 
beams. Thick beams composed of highly heterogeneous material lay-ups are considered. Comparisons are 
made with several beam theories, exact elasticity solutions, and results obtained with high-fidelity, two- 
dimensional elasticity finite element models. 


2. Basis for Zigzag Kinematics 


Laminated-composite beams are generally made of relatively thin layers of different materials, with 
different fiber orientations, that may have significantly different degrees of axial compliance. As a result, 
the axial displacement distribution of these heterogeneous, anisotropic beams exhibit varies in a nonlinear 
manner along the depth of the beam cross-section, especially for relatively deep beams with low 
slenderness ratios. Moreover, the axial displacement distribution is likely to exhibit a discontinuity in 
slope along the layer interfaces when the compliance of two contiguous layers is substantially different. 
Within each layer, the axial displacement distribution is presumed to be continuous, with continuous 
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slopes, consistent with the usual representation of the fibers and matrix material as a homogeneous layer. 
In a general layer-wise formulation of a beam theory, the axial displacement field could be modeled as a 
piecewise continuous function that is a collection of given nonlinear functions defined for each layer. 
However, if the beam stacking sequence is modeled as a collections of layers with layer thicknesses that 
are thin compared to the characteristic beam depth, then it is reasonable to presume that the nonlinear 
axial displacement distribution can be adequately modeled as a piecewise continuous function that is a 
collection linear functions defined for each layer. A function of this type is referred to herein as a 
“zigzag” function and, as a result, the corresponding beam theory is referred to as a Zigzag Beam Theory. 
Another important consideration, not addressed herein, is that if delaminations along given layer 
interfaces are present, a zigzag function with jump discontinuities at the given layer interfaces [20] could 
be used to represent the delaminations. 

Observations such as these have been substantiated by exact elasticity solutions such as that of Pagano 
[19]. Efforts like Pagano’ s prompted Di Sciuva [7] to formulate a refined plate theory in which a zigzag 
function is added to a special form of FSDT in which the shearing strains are used as independent 
variables. In this special formulation, the total rotations (the first derivatives of the deflection variable) 
provide the through-the-thickness linear distributions of the axial displacements, as in the classical 
Kirchhoff-Love plate and shell theories. Following Di Sciuva’ s work, Averill [14] proposed a similar 
zigzag formulation for the bending analysis of laminated-composite beams. In contrast, Averill’s 
formulation uses the standard form of TBT in which the bending rotation is the independent variable 
appearing in the kinematics for the axial displacement. However, both of these formulations have some 
undesirable attributes. In what follows, the essential aspects of the Di Sciuva and Averill zigzag models 
are examined in order to set the stage for the new refined zigzag beam-bending theory presented herein. 


2,1 Problem Formulation and Kinematics 


The basic problem of the present study is a narrow beam with depth of 2h and cross sectional area A 
that undergoes only plane deformations. The beam is made of N orthotropic-material layers that are 
perfectly bonded to each other, and the major principal material axis of each layer is parallel to the 
centrally located beam x axis shown in Figure 1. Static loads are considered which generally include a 
distributed transverse line load q(x) (force/length) and prescribed axial (T xa , T xh ) and shear (T za , T~ b ) 
tractions at the two reference cross sections, x = x a and x b , shown in Figure 1. Material points of the beam 
are located, prior to deformation, by the coordinates ( x,z ) , where x e \x a ,xf and z e [~h,h\ . 


A 



Figure 1. Beam subjected to transverse loading and end tractions. 
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For material points within the k -th layer of a laminated-composite beam, the through-the-depth 
coordinates are given by z e \z (k) ,z (k+]) \, where k = 1 . An expression for the beam displacements, 

that is general enough to describe the kinematics of the Di Sciuva, Averill, and present refined zigzag 
theories, is given by 

u (k \x,z ) = u{x) + z6(x) + </> {k \z ) y/(x) 
k z k) (x,z) = w(x) 

where the superscript ( k ) indicates *(*,^*(* + l) • In Eq. (1), the uniform axial displacement, u(x) , 
the bending rotation, 6(x) , and the transverse deflection, w(x) , are the primary kinematic variables of 
the underlying equivalent single-layer beam theory. The symbol (j) [k \z) denotes a piecewise linear 
zigzag function, yet to be established, and y/(x) is a primary kinematic variable that defines the 
amplitude of the zigzag function along a beam. When dynamic effects are considered, the four primary 
kinematic variables defined by Eq. (1) are also functions of time. If either (f) k ] (z)~ 0 or i//(x) = 0, the 
kinematic assumptions of Eq. (1) correspond to the TBT, if the rotation 6(x) is a primary kinematic 
variable. Note that depending on the specific FSDT theory used for the baseline kinematics, the four 
kinematic variables may have slightly different physical interpretations. 

The notations used in the present study together with a general distribution for the zigzag function are 
shown for a three-layer beam in Figure 2. The N + 1 axial displacements at the layer interfaces and the 
outer beam surfaces, as illustrated for TV = 3 in Figure 2, are denoted herein for an N -layer laminate by 
n {0 \u {1) ,...,ii {N) . Collectively, this set of displacements is referred to herein as the interfacial 
displacements even though u {0) and u [N) are not at the interface of two adjacent layers. 


1 

i 
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2h< 3 > 


k=3 


" 2h ,2> 


1c=2 

X 

l l 
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„ 2h< ’i 

l 
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(a) Layer notation for a three-layer laminate. 



(b) Zigzag function defined in terms of the interfacial 
axial displacements, u (t} ( i = 0 , 1 ,..., N) . 


Figure 2. Layer notation for a three-layer laminate and a corresponding generic zigzag function. 
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The linear strain-displacement relations of elasticity theory and the displacements defined by Eq. (1) 
give rise to the nonzero strain expressions 


4 k ) (x, z ) = u[ k l = u x + z 0 X + </> (k) y x 
y (k z \x, z) = u (k) : + = y + /3 {k) Y 


where /3 {k) = <f )' k) , and y = vv v + 0 represents the average shearing strain (or shearing angle) of the TBT. 

Note that, since (j> (k \z ) is piecewise linear, fi {k) is a piecewise constant function, i.e., it is constant 
across each material layer. For major principal material axes that are coincident with beam xaxis, 
Hooke’s stress-strain relations for the ^r-th orthotropic layer have the standard form 



where E (k) and G[ k) denote the axial and shear moduli of the k- th layer, respectively. 


To facilitate subsequent discussions, it is convenient to define a “difference function” rj(x) as 

J) = y~Y ( 4 ) 

This definition leads to the following convenient expressions for the transverse shear strain and shear 
stress in terms of the y(x) and rj(x) functions 

r«=(i +/?'*> )r (5) 

r «) = c m ( i +/? tt)) ) , _ G tty ‘> 7 ,6) 


2.2 Di Sciuva’s zigzag model 

Di Sciuva’s zigzag model [7] was originally developed for plate bending. The corresponding beam 
kinematics is obtained from Eqs. (1) by using the variable substitutions 

0 = y/-w x (7) 

and replacing z in Eqs. (1) with (z + h) . Thus, the axial displacement of Di Sciuva’s model has the 
form 


u (k) {x,z ) = u(x) + (z + h)[i//(x) - w x {x)] + (/) {k ^{z) y/{x) 


( 8 ) 



now 


where the transverse displacement is the same as in the second of Eqs. (1), and where <f>^ 

designates the actual form of the zigzag function used by Di Sciuva. These kinematics give rise to shear 
strain and stress that are piecewise constant (i.e., constant within each layer) and which are defined 
exclusively in terms of the amplitude function t// as follows 


A’=(i+/SV 


To determine the function (j )^ , Di Sciuva enforced a specific set of stress-continuity constraints along 
the layer interfaces that require all layers have the same (constant) transverse shear stress, that is. 


.(*) _ _(*+i) 

XZ L XZ ? 


k = 1, 


N-l 


( 10 ) 


Since Eqs. (10) impose only N — l constraints, and there are N + I interfacial displacements needed to 
define (f)^ , the zigzag function is specified to vanish across the entire bottom layer ( £/"" = u 1 1 =0), as 

illustrated in Figure 3(a). More generally, it is straightforward to select any given layer in which the 
zigzag function is specified to vanish [20]. Henceforth, determining a zigzag function by specifying (or 
fixing) a single layer contribution to vanish is referred to herein as a fixed-layer zigzag function method. 

The transverse shear stress distribution that results from Di Sciuva’s kinematics is uniform through the 
beam depth. If fixing the first layer ( k = 1 ) is used to define the zigzag function ( (f> ( r \ ] s = (3^1 = 0 ), as 
depicted in Figure 3(a), then the definition of t// as the shear strain in the first layer becomes evident 
from Eqs. (9); that is, 

w = rT (ii) 

In addition, the transverse shear stress in all layers simply equals that in the first layer and are given by 

C-<%r2 = (i2) 

An undesirable consequence of Di Sciuva’s approach to defining the zigzag function is the 
counterintuitive result that the transverse shear stresses in all layers are identical and equal to that in the 
layer that is specified to have zero- valued interfacial displacements. Moreover, the shear modulus of this 
fixed layer serves as a single weighting coefficient for the entire shear strain energy, which produces a 

bias in the deformation model toward the shear stiffness of the fixed layer. In contrast, (fi/fi should 
depend on all shear moduli and layer thicknesses, 2h {k) . The key property of Di Sciuva’s model is 

that the zigzag function vanishes identically for homogeneous beam, in which case the theory reverts to 
the underlying shear-deformation theory. The improved predictive capability of the model due to the 
</>ps \j/ term in Eq. (1) is particularly appreciable for thick and highly heterogeneous laminates. 


9 




(a) Di Sciuva’s zigzag-function distribution [7]. 



Figure 3. Di Sciuva’s and Averill’s zigzag-function distributions. 


Di Sciuva’s zigzag model exhibits another theoretical difficulty associated with the physical 
interpretation of the uniform transverse shear stress distribution that results from the procedure that 
defines its zigzag function. In particular, the physically correct shear force V x (x) obtained from the well- 
established equilibrium equation for beams 

V x = M XJt (13) 

is different from the corresponding force that is obtained when the shear stress in Eq. (12) is integrated 
over the cross-section, that is, 

V,*^'dA=Ar"' (14) 

A 

A related difficulty arises for a clamped support condition. For this case, the variationally consistent 
displacement boundary conditions are given as 

u = w = w x =i// = 0 (15) 

Because y/ = 0 is required at a clamped end, the corresponding transverse shearing strain and stress 
defined in Eqs. (9) vanish identically. In contrast, a non vanishing average shear stress T x , = V X / A at a 
clamped end is given by the shear force obtained from the equilibrium equation given by Eq. (13). 

From the perspective of finite element approximations appropriate for this type of beam theory, w(x) 
needs to be at least C 1 continuous, since the axial strain is proportional to the second derivative of w(x) . 
This requirement yields an added impediment for the approximation functions, particularly in developing 
efficient plate and shell elements based on this type of theory. It has been shown, however, that C° 
continuous kinematic approximations, usually associated with the FSDTs of Timoshe nk o [1], Mindlin 
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[17] and Reissner [18], result in simpler, computationally more efficient, and better performing finite 
elements than comparable elements based on C 1 continuous interpolations. It is this practical aspect that 
motivated Averill [14] to introduce a different zigzag model. 


2.3 Averill’s Zigzag Model 


Averill [14] developed a zigzag model for laminated-beam deformation that utilizes a penalty method 
that enforces the continuity of the transverse shear stress across the laminate depth in a limited sense. In 
this model, the standard form of TBT is used to provide the kinematics of the underlying equivalent 
single-layer theory; that is, the bending rotation is a primary kinematic variable in the expression for the 
axial displacement. The specific expression for the axial displacement used in Averill’s model may be 
expressed as 

u[ k \x,z ) = u{x) + {z + h-2h {l) )6{x) + </> A \z) y/( x ) (16) 

where (f)'\ 1 denotes Averill’s zigzag function, which is depicted in Figure 3(b) for a three-layer laminate. 
This expression is consistent with a zigzag function defined by fixing the interfacial displacements of the 
second layer ( k = 2 ). In Eq. (16), 0(x) represents the rotation of the second layer, that is, u[ 2) z = 6 . 

The shear stress continuity conditions given by Eq. (10) are enforced explicitly only on the first term 
of the shear stress given by Eq. (6). This enforcement is done by defining G A = G[ k, (\ + ff ! ’) and 
requiring this quantity to be constant for all material layers. Because the interfacial displacements of the 
second layer are fixed to define the zigzag function, it follows that (f ) 2 1 = fi {2 1 = 0 . Thus, it also follows 

that G a = G l2) , and the transverse shear stresses reduce to 


*™=G<?r -G^PVrj (k = \,...,N) 


(17) 


where y = y {2) , the shearing strain in the second (fixed) layer. Further, the second term in Eq. (17) 
diminishes to zero, in a limiting sense, by enforcing 

ri = y-y/^0 (18) 

with the aid of the following penalty constraint term which is added to the strain energy 

7 , f" V 2 dx (19) 


In this equation, A f is a penalty factor (units of force) which is set to a large value ( A a — » oo ) to ensure 

the validity of Eq. (18). Under these conditions, the shear stresses are uniform only in a limiting sense, 
that is, 


T 0c) Y+1) 

xz xz 


G™r (* = 1 , 


N- 1) 


( 20 ) 
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In addition, the transverse shear force derived from the beam equilibrium equation is not equal to the 
cross-sectional integral of the transverse shear stress, and is given as 

V x = jr^dA + k A 77 O) ( 21 ) 

A 

Averill’s zigzag model results in a four-variable theory and C°- continuous kinematics that are 
suitable for efficient finite element interpolations. However, the model appears to possess some 
significant theoretical anomalies. As in Di Sciuva’s zigzag model, there is the anomaly associated with 
specifying the interfacial displacements of a given layer to vanish, which is manifested as a bias toward 
the shear properties of that layer. This bias means that the methodology for defining the zigzag function 
lacks invariance with respect to the choice of the fixed layer. Moreover, Averill’s model gives erroneous 
results at a clamped end. More specifically, because a clamped boundary condition requires 


u = w = 0 = 1 // = 0 ( 22 ) 

and the penalty condition given by Eq. (18) yields the additional constraint w x — >0, the following 
erroneous conditions for the shear strain and stress at a clamped end result 

(/xz\ } )@clamped “ ^ 0 ( 22 ) 

support 

Averill recognized this anomaly and proposed to avoid prescribing boundary conditions on the 1 // 
variable, thus violating the variational consistency of the theory. 

The theoretical difficulties associated with Di Sciuva’s and AverilTs zigzag models motivated the 
development a variationally consistent zigzag model for laminated-composite beams that is devoid of 
these difficulties. The resulting new theory is presented subsequently and is referred to herein as the 
refined zigzag model. 


3. Refined Zigzag Model 


In this section, a novel, refined zigzag model for laminated-composite beams is presented that is free 
of the theoretical difficulties of the Di Sciuva and Averill models. The new model uses the TBT as the 
underlying equivalent single-layer theory that is obtained when the new zigzag function vanishes. The 
kinematic assumptions of the new model are given in Eqs. (1), and the corresponding strains and stresses 
are given in Eqs. (2)-(6). The key advantages of the new refined zigzag model are described subsequently, 
and the variationally consistent equilibrium equations and boundary conditions are derived from the 
virtual work principle. 


3.1 Zigzag Function 


The main difference between the new zigzag function presented herein and those used by Di Sciuva 
and Averill is that the new zigzag function is prescribed to vanish on the outer surfaces of the laminated 
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beam instead of within a given layer. This important difference ensures a more physically realistic 
distribution for the zigzag function. Specifically, the zigzag function distribution is nonzero throughout 
the beam depth and, as a result, the new zigzag function provides a contribution from each layer to the 
overall deformation of the beam cross-section. This new zigzag function is illustrated for a three-layer 
beam in Figure 4. 

With this new form of the zigzag function, the axial displacement given by the first of Eqs. (1) may be 
interpreted as a superposition of a layer-wise correction on to the standard form of TBT. That is, the 
terms u{x) and z 6(x) provide the basic, first-order, linear through-the-depth displacement distribution 
of TBT. Moreover, these terms provide the nonzero displacements at the bottom and top surfaces of the 
beam, which are zero-valued in the zigzag function; that is, 

u < x 0) = u(x) - h 6(x) and u[ N) = u(x) + h 6(x) (24) 

The zigzag component of the axial displacement, (f)'" ’i//(x) , represents the nonlinear distortion of the 
beam cross-section within the interior and is defined herein such that 

<f) ik) = (j)(z,u (k ’) (k = 1,..., N — l) (25) 

This definition for the zigzag function is intended to indicate that its through-the-depth variation depends 
only on the interior displacements at the layer interfaces, as illustrated in Figure 4 for a three-layer 
laminate. 



A 

r 



■' 2H 3) 


k=3 


" 2h< 2 > 


1c=2 

X 

l l 



W 

, 2hP\ 


k=l 



(a) Layer notation for a three-layer laminate. 



(b) Zigzag function of the refined zigzag theory 
defined in terms of the interfacial axial displacements, 
u u) (i = 0,l,..., N). 


Figure 4. Layer notation for a three-layer laminate and zigzag function of the refined zigzag theory. 


With this type of kinematic arrangement, the axial strains at the top and bottom beam surfaces are 
represented by the TBT contribution to the axial displacement field; that is, 


13 



(26) 


= U x - h 6 x and s[ N) = u x + h 6 x 

To arrive at a simple representation of the zigzag function, (f) {k) , in terms of the layer-interface 

displacements, u (k 1 , it is convenient to introduce local layer coordinates, tf kj , that are non- 
dimensionalized by the corresponding layer thickness, 2h {k) , such that % {k) e [-1,1] and k = l,...,N. 
These local thickness coordinates are defined in terms of the N -layer-laminate thickness coordinate 
z e \-h,h\ as 


g( , ) _ z-(z^ + ^>) 


( k = \,...,N ) 


(27) 


where z <0) = —h, z {k> = z (k ” + 2h ik ’ , and z <N> = h. For a linear variation of the zigzag function across 

the thickness of every layer, with ($ k> = i/ k l> and ($ k) = u (k> for ^ k> = - 1 and ^ k> = 1 , respectively, the 
zigzag function used in Eqs. (1) is given by 




(28) 


where u i0) = w <,V) = 0. The above form of <j) (k) gives rise to a piecewise constant function fi (k) 
( p { k ] = (j) k 1 ) that may be expressed in terms of the layer-interface displacements as 


J3 (k) = 


u 


^-u (k - l) 

2h (k) 


(29) 


With this expression and the prescribed conditions u {0) = u iN) = 0 , it is seen that the function fT k ’ has a 
particularly useful property that is true regardless of the layer material properties and thicknesses; that is, 


f [3 {k) dA = — Y 2K lk) j8 { [k) = — U [N) - u m ) = 0 

j 2 h k % 2h { ’ 


(30) 


The usefulness of Eq. (30) is seen by integrating Eq. (5) over the beam cross-section. Performing this 
process and then recognizing the above J3 {k ) property, gives 



This equation verifies that y(x) is in fact the shear angle of the TBT and represents an average shear 
strain of the cross-section. Another useful application of Eq. (30) is found by multiplying both sides of 
Eq. (5) by /3 {k) , and then integrating over the cross-sectional area. This process reveals that the function 
y/(x) is a weighted-average shear strain quantity given by 
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1 


(32) 


' JO 6 w fdA 

A 


ffVJ’JA 


At this stage in the analysis, the layer-interface displacement, u (k) , are yet u nk nown. Recall that in the 
zigzag models of Di Sciuva and Averill, the zigzag function was defined by prescribing the layer- 
interface displacements to be zero-valued for a specified layer, and by enforcing continuity of the 
transverse shear stresses at each layer interface. In the present study, instead of prescribing the layer- 
interface displacements to be zero-valued for a specified layer, the conditions u {0) = u ' = 0 have been 
used. This approach allows a physically attractive partitioning of the TBT and zigzag kinematics, without 
introducing a bias in the kinematic equations. 

In the manner of both the Di Sciuva and Averill schemes, the new refined zigzag model enforces 
interfacial continuity of the first term, associated with the /( X ) variable, in the expression for the 
transverse shear stresses given by Eq. (6). Applying this process yields the conditions 

G W(i +A W) =G (^ 1) (i + ^ + i) ) (33) 

for k = 1,..., N - 1 . Thus, Eqs. (33) effectively establish a beam material constant defined by 

G = (l + /3 {k) ) (constant) (34a) 

Whereas this constraint strategy is fully consistent with the shear-stress continuity constraints developed 
by Di Sciuva, Eqs. (9) and (10), in the present setting because of an additional kinematic variable ij , the 
conditions (33) do not enforce the interface shear-stress continuity. Furthermore, in contrast to Averill’s 
penalty formulation, which enforces 77 — >0(Eq. (18)) and thus achieves shear-stress continuity in a 
limiting sense, the present refined zigzag model does not seek nor enforce continuity of the transverse 
shear stresses along the layer interfaces. 

Obtaining /3 {k) from Eq. (34a), 

J3 (k) = G/G% ) -\ (34b) 


and then invoking Eq. (30), i.e., integrating Eq. (34b) across the beam cross-section, gives rise to the 
expression for G in terms of the shear modulus and thickness of every layer. The resulting expression is 


G = 




dA 

G^ , 

xz J 




z 


h (k) 


G=in v xz j 


(35) 


Then, a recursion relation for the layer-interface displacements is found by substituting Eq. (34b) into Eq. 
(29) and accounting for Eq. (35); that is, 
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u {k) = u (k 0 + 2h 


( k ) 


c (k) y 

xz 


Ak) 


-1 


% G 


(A 

XZ J 


; (u (0) = u (N) = 0; £ = !,.. .,7V) 


(36) 


By substituting Eq. (34b) into Eq. (5), and using Eq. (4), the transverse shear strains may be conveniently 
expressed in terms of the constant and piecewise-constant components according to 


7^ = V + 


;(*) 


¥ 


(37) 


Similarly, invoking Hooke’s relations, Eqs. (3), the transverse shear stresses may be expressed as 

A = G V + G«> n = G y + (G® - G) 7 (38) 

An alternative procedure that determines the same 0 (/l) function as the one just presented can be 
developed from the energy minimization point of view and is summarized in Appendix C. 


3.2 Constant stress requirement and penalization issues 


The key issues concerning the validity and theoretical difficulties of the Di Sciuva and Averill zigzag 
models can now be clarified from the insight provided by the new, refined zigzag model. First, recall that 
in the Di Sciuva [7] and Averill [14] models, enforcing continuity of the transverse shear stresses across a 
cross-section, whether explicitly or as a limiting condition, results in an erroneous constant value for the 
transverse shear stress in a layer that does not depend on the corresponding shear moduli. This 
representation of the transverse shear stresses is generally a poor approximation of the actual transverse 
shear stresses present in a beam with a heterogeneous cross-section. Moreover, their approach produced a 
physical anomaly for a clamped boundary condition. The new zigzag model and the results of the present 
study reveal that the anomalies associated with the Di Sciuva and Averill models arise because the stress 
continuity constraints are fundamentally inconsistent with the kinematic freedom inherent in the lower- 
order kinematic approximations of the underlying beam theory. Naturally, higher-order approximations, 
such as the cubic zigzag model in [ 12 ], may provide the kinematic freedom required for consistent, 
continuous, non-constant transverse shear stresses. 

The second issue that was clarified during the present study is the penalty enforcement, in AverilTs 
model, which also produced a theoretical anomaly for a clamped boundary condition. The difficulty with 
this particular boundary condition was resolved unsatisfactorily in [14]. During the course of the present 
study, it was found that the penalization issue associated with the f] function can be properly resolved at 
the constitutive-equation level, as opposed to adding an ad hoc penalty functional to the variational 
principle, as in [14]. In particular, a dimensionless penalty factor ,A 0 , can be introduced into the 
constitutive equation for the transverse shear stresses as follows 

=Gr + M<%’-G)i7 ( 39 ) 
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where, for the special case /t n = 1 , Eq. (39) reverts back to Eq. (38), which is the expression for the 
transverse shear stresses used in the refined zigzag model (although, for demonstration purposes, a few 
numerical results have been carried out with a large value of X 0 and demonstrated in Section 5). The 
transverse shear strain energy associated with Eq. (39) is given by 

u$r =\) Ki, y™ dAdx = ) y2 dx + f W dx (40) 

x a A x a x a 


where X denotes a penalty parameter (units of force) defined as 



dA-GA 

J 


(41) 


The first term in the shear strain energy has the same exact form as the shear strain energy of TBT, for a 
homogeneous beam, provided G is the beam’s shear modulus. The second term is associated with the 
penalty parameter, X . This second term appears naturally as a consequence of the modified constitutive 
relations and is amenable to physical interpretation. For homogeneous cross-sections, X vanishes because 
J3 ik) =0 and G reduces to the beam’s shear modulus. This approach is significantly different from 
Averill’s penalization approach, where a similar penalty term is added to the strain energy a posteriori. 
The inclusion of the penalty factor X 0 in Eq. (41), under the condition^ — » oo , allows for the limiting 
condition in Averill’s model, given by Eq. (18), to be fulfilled. Also, note that the coupling term yrj 
vanishes from the shear energy because of the intrinsic property of the new, refined zigzag function given 
by Eq. (30). 

The fully consistent form of the present zigzag model is given by X 0 = 1 in which Averill’s limiting 

condition, Eq. (18), is not enforced. As a result, using X 0 = l yields transverse shear stresses and strains 

that are piecewise constant across the laminated-beam cross-section and no anomalies are present for a 
clamped boundary condition. In addition, the same basic conclusions may be reached for any arbitrary 

value ofX 0 , including very large values. In this case, // = 0( A 1 ) , where O denotes the order of the 

quantity in parenthesis. Even for values ofX 0 -A oo , the second term of the shear stress expression in Eq. 

(39) does not vanish and is 0(G- G [ k) ) . 


3.3 Virtual work principle, equilibrium equations, and boundary conditions 

For the beam defined in Figure 1, the virtual work principle (VWP) may be expressed as 
J J(°f )& i* : ’ + T^Sy^dAdx- \qSwdx 

* a A X a (42) 

+ ^T m 8u*\x a ,z) + T m 8vix a j\dA- \\T xh Su { x k) (x h ,z) + T b Sw(x b )~\ dA = 0 


17 



where 8 denotes the variational operator. Performing the cross-sectional integration and variation by 
parts, results in the following equilibrium equations 


K, x ~v x = o 
v,,+q = o 
m^-v 4 = 0 


and a consistent set of boundary conditions 


either u(x a ) = u a 
either 0(x a ) = 0 a 
either w(x a ) = w a 
either V /(x a ) = y7 a 


or N x (x a )=N xa 
or M x (x a )= M xa 

or K(xJ= K a 

or M^x a )=M^ a 


(43) 


(44) 


where a=a,b and the bar-superscripted symbols denote the prescribed displacements and stress resultants. 


The first three equilibrium equations in Eqs. (43) are the same as those in the conventional Euler- 
Bemoulli Beam Theory and Timoshe nk o Beam Theory. The fourth equation describes the moment-shear 
equilibrium associated with the contribution of the zigzag function to the bending moment and shear 
force. The stress resultants appearing in these equations are given by 




K 


At, 


K, ^ ^)dA 


A 

(v„, M„, A/ # , F„)=J(r„, zT„, %, T ;a )dA 

A 


(45) 

(46) 


where N x , M x , and V x are respectively the conventional axial force, bending moment, and shear force. 

The symbols M , and V, are the bending moment and shear force associated with the contribution of the 

zigzag function. The corresponding beam constitutive relations are obtained by substituting the strain- 
displacement relations into the constitutive relations, given by Eqs. (3), and then substituting the resulting 
expression into the stress resultants defined by Eqs. (45) and integrating over the cross-sectional area. The 
resulting constitutive relations of the new zigzag beam theory are 


K 


A\ Ai 


U , 

K 

>= 

B\2 E\\ ^2 

£ 

M, 



Y*. 



of 

\ 0+W , x 


Qn Q 22 


V 


(47) 


(48) 
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where the stiffness coefficients are defined in Appendix A. 


Substituting the expressions for the stress resultants in Eqs. (47) and (48) into Eqs. (43) yields the 
equilibrium equations in terms of the primary kinematic variables of the theory. This eighth-order system 
of ordinary differential equations, with constant coefficients, is given by 

A AlAxt AW.xx — ^ 

Ql: \ ( w ,.u + #v) +Qk.v =~^ x ) 

B \2 U xx + D lAxx + D n'/ / xx -Ql ( W x +^) -Ql Y = ® 

B \i u ,xx + D vAxx + D n¥,xK - Ql ( w ’r +#)“ QiV = 0 

Solution of the corresponding boundary-value problem involves integration of the four equilibrium 
equations, given by Eqs. (49), subject to the boundary conditions given by Eqs. (44). For problems with 
forces specified at the ends of a beam, the constitutive relations given by Eqs. (47) and (48) are used to 
express the force boundary conditions in terms of the primary kinematic variables. 


4. Examples and Analytic Solutions 


The predictive capability of the new, refined zigzag beam theory is demonstrated for two beam 
examples that provide a challenge to the theory for certain laminate constructions and beam-slendemess 
ratios. The first example is a simply supported beam of length L that is subjected to the sinusoidal 
transverse line load q(x) = q 0 sin ( nx / L ) , as shown in Figure 5. 



Figure 5. Simply supported beam subjected to a sinusoidal transverse line load. 


The boundary conditions at the ends x = 0 and x = L of the beam are 
w = N x = M x = M^ = 0 


( 50 ) 


The following trigonometric expressions for the four kinematic variables 
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w(x) = w 0 sin(;rx/Z), 

( u{x ), <9(x), j/(x)) = (w 0 , 0 O , ^ O )cos(;rx/Z) 


( 51 ) 


satisfy the above boundary conditions exactly. The amplitudes w 0 , u 0 , 0 O , and t// 0 are determined 

uniquely, and the boundary-value problem is solved, by substituting Eqs. (51) into the equilibrium 
equations, Eqs. (49), and then solving the resulting system of four linear algebraic equations. 

The second example presented herein is a cantilevered beam of length L that is subjected to a 
transverse load F at the free end, as shown in Figure 6. The boundary conditions are given by 

x = 0: u : = w= 6 = y/ = 0 

(52) 

x = L : N x =M x =M^ = 0, V x = F 

Because q(x ) = 0 for this problem, the equilibrium equations, Eqs. (49), can be rewritten in the more 
convenient form 


£>, x,, + Q//., = 0 
c 

^xx=-Q'/ / xv+7jf'/ / 

Ml 

w , x =-0- C a¥ + C 5 d xx + c b \f/ xx 

u =—C^0 —C„y/ 

,xx 7 ,xx 8 T ,xx 


(53) 




L, 






Figure 6. Cantilevered beam subjected to a transverse tip load. 


The coefficients, which are dependent on the material properties, are presented in Appendix B. 
Integrating these equations by using the standard technique for ordinary differential equations yields the 
following expressions for the kinematic variables 

y/{x) = a, cosh(itr) + a 2 sinh(itr) + a 3 , 
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0(x) = 


( c ^ 

-c 3 + -fV 

R 2 D nJ 


(a t cosh(i?x) + a 2 sinh(itt)) + 


Cm 


7 W 1 9 

x + a 4 x + a 5 , 


2D 
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u{x) = 


-c, + c,c, - 


cc. 
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2 '-'7 


8 8 7 „2 n * 




C* C a 

(a, cosh(itc) + <x, sinh(i?x)) 2 7 , 3 x 2 + a 6 x + a 7 , (54) 

2D n 


w(x) = 


Q c i 
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C 2 C 5 r 

d: a 


+r(c 6 -c 3 c 5 ) 


(«j sinh(i?x) + a 2 cosh(i?x)) 


C~,a 3 1 rt. 2 
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6D U 
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v u \\ 


a 3 —a 5 


x + a„ 


where the a i (i = 1,...,8) u nk nown constants are determined from the boundary conditions stated by Eqs. 
(52). 


5. Results and Discussion 


Results are presented in this section for the two beam examples previously described. The results 
include comparisons of those obtained from the new zigzag theory to those obtained from the TBT, Di 
Sciuva’s and AverilTs zigzag models, exact elasticity solutions, and high-fidelity finite element solutions. 
A thick, three-layer beam is considered for all problems, with total depth 2h = 2 [cm] and length L = 10 
[cm]. The span-to-depth ratio is L/2h = 5 and each beam has a rectangular cross-section. The 
mechanical properties of six layer materials that were used to generate results are presented in Table 1. 
The major principal axis of each material is aligned parallel to the beam axis. Three unidirectional 
laminate stacking sequences that were considered are presented in Table 2. 


Table 1. Mechanical properties of the beam layer materials. 


Material 

E[ k) [GPa] 

Git’ [GPa] 

type 

(Young’s modulus) 

(Shear modulus) 

a 

73.0 

29.2 

b 

21.9 

8.76 

c 

3.65 

1.46 

d 

0.73 

0.292 

e 

0.219 

0.088 

f 

0.073 

0.029 
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Table 2. Stacking sequences of three-layer beams 
(layer sequence is in the positive z direction). 


Laminate 

Thicknesses [cm] 

Materials 

A 

(0.20/1.60/0.20) 

(a/f/b) 

B 

(0.66/0.66/0.66) 

(b/e/b) 

C 

(0.66/0.66/0.66) 

(d/a/c) 


The layer thicknesses are presented in the form {2h (X) , 2ti’ ] , 2h <} ) and the first layer in each laminate 
starts at z = —h. Similarly, the material composition of each layer is presented in the form 
( /V/"\ M {2 \ M 1 ! 1 J , where “M” denotes the material type. The first stacking sequence, designated as 

Laminate A in Table 2, is an asymmetrically laminated sandwich construction with a core that is eight 
times thicker than the face sheets. In addition, the core is made of material “f “ and is three orders of 
magnitude more compliant than the bottom face sheet, which is made of material “a.” Moreover, the top 
face sheet has the same thickness as the bottom face sheet, but is about three times more compliant. The 
extreme heterogeneity of this stacking sequence results in beam examples that are very challenging for 
the beam theories considered herein to model adequately. Laminate B has three equal-thickness layers, 
with the outer layers having the same stiffness properties, and where the middle layer is two orders of 
magnitude more compliant than the outer layers. Moreover, this laminate is symmetric with respect to the 
mid-depth reference axis. Laminate C also has three equal-thickness layers; however, unlike Laminate B, 
here the middle layer is two orders of magnitude stiffer compared to the bottom layer. Also note that this 
laminate does not possess material symmetry with respect to the mid-depth reference axis since the top 
layer, made of material “c”, is five times stiffer than the bottom layer, made of material “d”. These three 
distinctly different laminates thus constitute a diverse set of practical and challenging material 
configurations that represent realistic bounds in industrial applications. 

Results for simply supported beams made from Laminates A and C are presented in Figures 7-12, and 
results for cantilevered beams made from Laminates A and B are presented in Figures 13-26. Several 
curves are presented in each figure, as indicated by their figure legends. For conciseness, the following 
legend captions are used. The legend caption “Exact” refers to the exact elasticity solution for a simply 
supported plate in a state of cylindrical bending that was presented by Pagano [19]. Although this solution 
is for a plate, the results are also applicable to beams. The legend caption “FEM/NASTRAN” refers to a 
high-fidelity, two-dimensional FEM solution obtained with the MSC/NASTRAN® code [21]. The details 
of this solution are described in Appendix D. The legend caption “TBT” refers to the exact solutions 
obtained from the Timoshenko Beam Theory. For these solutions, the standard shear correction factor 

k 2 =5/6 was used throughout, which is appropriate for homogeneous beam with a rectangular cross- 
section. Note that various approaches have been proposed for determining shear correction factors for 
laminated beams that often provide more accurate transverse deflection results. The choice of a shear 
correction to be used with TBT, however, does not influence the prediction of the axial and transverse 
shear stresses. The legend caption “Zigzag (D, A)” refers to exact solutions obtained with the Di Sciuva 
[7] and Averill [14] zigzag models, for the cases when virtually identical results are produced. The 
caption “Zigzag (R)” indicates solutions obtained using the refined zigzag theory with Z 0 =l. Similarly, 
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“Zigzag (R(/t w =l O 6 ))” refers to solutions obtained using the refined zigzag theory with Z 0 =10 6 . The 
caption “Zigzag (All)” refers to nearly identical solutions obtained by all three zigzag models presented 
herein. Lastly, the legend caption “Integrated” refers to the transverse shear stresses obtained by 

integrating the equilibrium equation a[ k ' + T [k _\ = 0 , in which u (k) represents the axial stress determined 
from the refined zigzag theory (A#= 1). This commonly employed transverse shear stress-recovery method 
has been shown to produce accurate shear stresses T ( x k> that are in close agreement with the corresponding 
results obtained from the high-fidelity finite element model described in Appendix D. 

In Figures 7-26, the axial and transverse coordinates are normalized as (x, z) = (x/ L, z/h). 
Similarly, for the simply supported beams, the following dimensionless solution variables are used 




10 q 0 L 


m L‘>) 


5 ^ XI 


q 0 L 


For the cantilevered beams, the dimensionless solution variables are given by the expressions 

(s'T. 

K=~\ t T‘U 

x p Ja xz 


ul k) , w 


\ - Ai (..(k) <k) \ (~(k) ~(t)\ 2hA / (t) 

' = ToaA ' x ’ z ’ ‘ r ’ - !~fT' x ’ xz ' 


(55) 


(56) 


5.1 Results for Simply Supported Beams 


Normalized axial displacements, transverse displacements, and transverse shear stress distributions are 
presented in Figures 7-9, respectively, for the simply supported beam made from Laminate A - a three- 
layer sandwich with a very thick, compliant core. Three curves are shown in Figures 7 and 9 and two 
curves are shown in Figure 8. The solid and dashed black curves correspond to results from the exact 
elasticity solution and all three zigzag models presented herein, respectively. The finely dashed gray line 
corresponds to results from the TBT. The abscissa in Figure 8 is scaled to exaggerate the difference in the 
results obtained from the exact elasticity solution and the three zigzag models. As a result, the transverse 
deflection obtained from the TBT, u z (L/2) = 0.164 , is not shown. 

The results in Figures 7-9 show that modeling the deformations of Laminate A is a particularly 
challenging problem for any structural theory. The beam is relatively deep (L / 2h = 5) and has a high 
degree of heterogeneity typical of sandwich construction: a highly compliant thick layer bounded by two 
relatively stiff, thin face sheets. For this problem, Timoshe nk o theory under predicts the deflection by a 
factor of about 15 (Figure 8). All zigzag theories produce comparable displacement and stress results that 
correlate well with the exact solutions. As seen from Figure 9, only slight differences are observed in the 
shear stress as predicted by the Di Sciuva, Averil and refined zigzag theories, its distribution being 
respectively constant and piece-wise constant. Contrasted with the zigzag results, Timoshe nk o’s shear 
stress is substantially under estimated in the core and grossly over estimated in the face sheets. 

Laminate C has three layers of equal thickness, with a very stiff inner layer (Figures 10-12). For this 
laminate, the TBT models the overall beam stiffness adequately and produces accurate displacement and 
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transverse shear stress predictions. As shown in Figure 11, the refined zigzag theory, when used with a 
large penalty factor, A f) = 1 0 6 , produces a slightly lower deflection than in the A n = 1 case, where the latter 
solution only slightly exceeds the exact value atz = 0. Also note that the deflection solution 
corresponding to the Ao= 1 (/’ case is closely correlated with the Di Sciuva and Averill predictions that are 
somewhat stiffer than those of the refined theory with Ao= 1 . As evidenced from Figure 12, the refined 
zigzag theory demonstrates a notable improvement in the transverse shear stress predictions over the other 
zigzag theories considered herein. This example clearly demonstrates that, even for the simply supported 
case, the other zigzag theories may produce rather inaccurate shear stresses, even when compared to those 
of the TBT. Moreover, the non penalized form (T 0 =l) of the refined zigzag theory produces slightly more 
flexible and accurate results than its penalized (A 0 = 10 6 ) counterpart. For this reason, subsequent results 
for the refined zigzag theory are only demonstrated for the non penalized case (Ao= 1). 



Figure 7. Normalized axial displacement at the left end of a simply supported beam with 

laminate A stacking sequence. 
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Figure 8. Normalized transverse mid-span displacement of simply supported beam with 

laminate A stacking sequence. 
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Figure 9. Normalized transverse shear stress at left end of simply supported beam with 

laminate A stacking sequence. 






Figure 10. Normalized axial displacement at left end of simply supported beam with laminate C 

stacking sequence. 



Figure 11. Normalized transverse displacement at mid-span of simply supported beam with 

laminate C stacking sequence. 
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Figure 12. Normalized transverse shear stress at left end of simply supported beam with laminate 

C stacking sequence. 


5.2 Results for Cantilevered Beams 


In this section, results are examined for the cantilevered beam problem for Laminates A and B 
(Figures 13-26). The sandwich-type Laminate A is again a major challenge for the TBT, which over 
estimates the over all stiffness of the beam and under estimates the maximum deflection by over 80% and 
the maximum axial stress by over 50% (Figures 13 and 14). These results are in great contrast with those 
of the zigzag theories, which produce very accurate displacements and stresses. 

Differences between the various zigzag solutions are observed by examining the through-the-depth 
distributions predicted for the displacements and stresses at various locations along the span, including 
the clamped end. The axial stress and displacement depicted in Figures 14-16 are accurately modeled by 
all three zigzag theories examined, exhibiting only minor quantitative differences. The good agreement is 
attributed to their built-in zigzag kinematics, giving rise to the requisite slope changes at the layer 
interfaces due to the layer differences in shear moduli. 

Examination of the transverse shear stress predictions reveals specific differences in the capabilities of 
the beam theories examined. At the clamped end (Figures 17), the refined zigzag theory yields a non- 
vanishing stress, whereas the Di Sciuva and Averill zigzag theories predict an erroneous, zero shear 
stress. Away from the clamped end (Figures 18 and 19), the shear stress distributions also exhibit 
substantial differences. For this problem, the cross-sectional distribution of the shear stress changes 
significantly along the beam's span, as expressed by the shear stress obtained through the integration 
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procedure (highlighted in Section 5) and the refined zigzag theory. The Di Sciuva and Averill zigzag 
theories produce a shear stress that varies in magnitude with respect to the axial coordinate, while 
remaining constant across individual cross-sections. The shear stress of the refined zigzag theory also 
varies with respect to the axial coordinate; however, it is piece-wise constant across the cross-sections. In 
contrast, the shear stress of the TBT remains constant with respect to the axial coordinate and exhibits 
inferior accuracy across the beam’s cross-sections. Finally, the shear force distribution along the beam’s 
span (Figure 20) clearly illustrates the transverse-shear force anomaly inherent in the Di Sciuva and 
Averill zigzag theories. As shown in Figure 20, the Di Sciuva and Averill zigzag theories predict that the 

normalized shear force quantity, V x , vanishes at the clamped end. Nearly halfway across the beam, an 

asymptotic value is achieved that over estimates the correct value by about 10%. Note, however, that a 
correct shear force distribution is obtained from Di Sciuva’s theory when it is calculated by taking the 
derivative of the bending moment. Both the TBT and the refined zigzag theory give rise to a correct, 
constant shear force over the entire span of the beam. 

Laminate B has three equal-thickness layers and possesses a material symmetry with respect to the 
mid-depth reference axis. For this laminate, the solution obtained by the TBT is overly stiff, with the 
maximum deflection underestimated by about 50% (Figure 21). In contrast, results from the refined 
zigzag theory are consistently accurate. Specifically, the refined zigzag theory models the span-wise 
variation in the transverse shear stress distribution across each cross-section accurately (Figures 22-25). In 

Figure 26 are depicted the normalized shear force distributions, V x , along the beam span. As in Figure 
20, the key theoretical anomaly of Di Sciuva’s theory is manifested by V x vanishing at the clamped end. 

It is noted once again that a correct shear force distribution corresponding to Di Sciuva’s theory is given 
by the derivative of the bending moment. 



Figure 13. Normalized deflection of cantilevered beam with laminate A stacking 

sequence. 
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Figure 14. Normalized axial stress at the clamped end of cantilevered beam with laminate 

A stacking sequence. 



Figure 15. Normalized axial displacement at the free end of cantilevered beam with 
laminate A stacking sequence. 






Figure 16. Normalized axial displacement along the top face of a cantilevered beam with 

laminate A stacking sequence. 




Figure 17. Normalized transverse shear stress at the clamped end of a cantilevered beam with 

laminate A stacking sequence. 


30 







Figure 18. Normalized transverse shear stress at x=L/10 of a cantilevered beam with laminate 

A stacking sequence. 



Figure 19. Normalized transverse shear stress at the free end of a cantilevered beam with 

laminate A stacking sequence. 
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Figure 21. Normalized deflection along a cantilevered beam with laminate B stacking sequence. 








Figure 22. Normalized transverse shear stress at the clamped end of a cantilevered beam with 

laminate B stacking sequence. 




Figure 23. Normalized transverse shear stress at x=L/10 of a cantilevered beam with laminate B 

stacking sequence. 
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Figure 24. Normalized transverse shear stress at x=2L/5 of a cantilevered beam with laminate 

B stacking sequence. 




Figure 25. Normalized transverse shear stress at the free end of a cantilevered beam with 

laminate B stacking sequence. 
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Figure 26. Normalized shear force along the span of a cantilevered beam with laminate B 

stacking sequence. 


6. Conclusions 


A new, refined “zigzag” theory for laminated-composite and sandwich beams that exhibit a high 
degree of transverse shear flexibility has been presented. The kinematic assumptions of this theory 
employ a novel piecewise continuous representation of the displacement field that is independent of the 
number of material layers, does not require enforcement of transverse-shear-stress continuity to yield 
accurate results, and that is zero-valued at the top and bottom surfaces of a beam. Additionally, the force 
equilibrium equations, boundary conditions, and strain-displacement relations are completely consistent, 
with respect to the virtual work principle, and a transverse-shear correction factor is not required. This 
new theory is better suited for engineering practice than previous theories because of its relative 
simplicity and its ability to model accurately the transverse shear and axial deformations of the individual 
layers in a physically realistic manner. Unlike other theories, meaningful axial and shear stresses are 
obtained directly from the constitutive equations, in a theoretically consistent manner, in this new theory. 
Moreover, the resultant transverse shear force acting at a beam cross-section can be obtained by 
integrating the corresponding transverse shear stress directly. These two physically attractive, 
fundamental properties are unique to this new theory. Similarly, this new theory is devoid of a major 
shortcoming of other similar theories; that is, the new theory enables accurate modeling of the clamped 
boundary condition. 
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Results have been presented for simply supported and clamped beams with highly heterogeneous 
material lay-ups and relatively small span-to-depth ratios. These results are for beam constructions that 
are extremely difficult for beam theories to model accurately. Comparisons with corresponding results 
obtained from exact two-dimensional elasticity analyses, high-fidelity finite element analyses, and other 
beam theories have also been presented. These comparisons show that the new, refined zigzag theory 
provides highly accurate predictions of displacements, stresses, and force resultants. 

Finally, this new theory is amenable to the development of computationally efficient (^-continuous 
finite elements that can be implemented into general purpose codes used to design an important class of 
high-performance aerospace structures, and is readily extendable to theories for thick, highly 
heterogeneous plate and shell structures. 
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Appendix A. Stiffness coefficients in refined zigzag theory 


Listed below are definitions of the stiffness coefficients that appear in the constitutive equations, Eqs. 
(47) and (48). 


C,„-b, 2 .a,)=£ Ey(\,zy)dA 

(Al) 

(B a ,D n ,D 21 ) = lEyt m (\,z,t m )dA 

(A2) 

(afe.feiHe+A-'M) 

(A3) 

Q = [ G (k) (1 + p (k) f dA = GA 

JA 

(A4) 


where the constants G and X are defined in Eqs. (35) and (41), respectively. Note that a Timoshenko- 
type shear correction factor is not used within the present theory when modeling a heterogeneous cross- 
section. For homogeneous beams or laminated beams in which the shear moduli are the same (or nearly 
the same), the present theory reduces back to Timoshe nk o theory. In this case, it is appropriate to use the 

standard shear correction factor k 2 =5/6 (for a rectangular cross-section) in the definition of the shear 
stress. 
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Appendix B. Cantilevered beam: Exact solution material coefficients 


The exact solution in Eqs. (54) corresponds to the case C 2 /Cj< 0. For C/C;>0, trigonometric functions 
define the exact solution (instead of hyperbolic functions), in which case some coefficients listed below 
would have different expressions. 
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Appendix C. Alternative derivation of zigzag function 


An identical form of the zigzag function </> {k) derived in Section 3.1 can be obtained by minimizing 

the shear strain energy for the case 77 = 0. Under this assumption, the transverse shear strain energy per 
unit length and width may be written as 


•s = \ ] G2> (1 + P'")r (1 + P m )r dz = (rF G " > ( 1+ / ?m ) ! lh 

2 _ h 2 k = 1 
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Using Eq. (29) one obtains 
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Minimizing S with respect to the u nk nown interface displacements u (k> yields 
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resulting in (for k=l , ..., A-7) 
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(C4) 
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The above relationship is identical to that of Eq. (33). 


Appendix D. Reference finite element model for the cantilevered beam 


The cantilevered beam is analyzed using a two-dimensional plane-stress model discretized with 
QUAD8 elements based on serendipity shape functions [21]. The discretization involves a uniform mesh 
having 40 elements in the direction of the beam thickness and 200 along its length, with a total of 8000 
elements, and nearly 49,000 degrees-of-freedom. The tip transverse force F is modeled by way of a 
parabolic shear traction f(z ) ( / (+/?) = 0 ) distributed along the free end of the beam, producing the 

total shear force of F = f(z)dz . The displacements are computed at the nodes, whereas the strains 

J-h 

and stresses are evaluated at the element centroids. 
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Figure 27. Geometry, boundary conditions, and loading for the cantilevered beam modeled as a two-dimensional 

plane-stress problem using MSC/NASTRAN® 1 . 
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